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Abstract 


Operational remote sensing of terrestrial atmosphere is heavily based on the 1-D radiative transfer 
equation. However, cloudy scenes are influenced by 3-D effects (e.g., illumination from cloud sides, etc.). 
This leads to biases in aerosol/cloud/land/ocean retrieval schemes for scenes with clouds. These biases can 
be understood and quantified only with the use of the 3-D radiative transfer theory, which allows to 
account for arbitrary spatial variation of atmospheric parameters. The task of this paper is twofold. First of 
all we introduce a novel technique for the solution of the 3-D radiative transfer equation based on the grid 
approximations and the straightforward iteration procedure realised on supercomputers with parallel 
architecture. We study the performance of our technique comparing with the solutions obtained by the 
Monte-Carlo code. A close correspondence is found. Secondly, we quantify the influence of neighbouring 
clouds on the clear sky reflection function at the nadir observation depending on the solar illumination 
conditions. We find that the influence of cloud on the clear sky reflectance function is not negligible (even 
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outside the cloud geometrical shadow). Thus, the peculiar inner boundary layer arises in the sky reflectance 
function with shadowing and brightening effects. 


© 2004 Elsevier Ltd. All rights reserved. 
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1. Introduction 


Radiative transfer in a cloudy atmosphere is usually studied in the framework of the plane- 
parallel approximation. Then the diffused light field changes only along the vertical direction for a 
wide solar beam illumination conditions. There is no change in the radiation field in the horizontal 
direction. Although this approximation, which is often called 1-D case, is very important for the 
radiative transfer studies of extended cloudiness (e.g., extended fields of Stratocumulus clouds), it 
cannot be applied for majority of cloudy scenes. Indeed cloudiness has a horizontal structure (e.g., 
holes between clouds). 

These effects can be accounted for in the framework of the 3-D radiative transfer 
equation, where the spatial variation of local optical properties is fully accounted for. Various 
approaches to deal with 3-D clouds are known [1]. Most popular techniques are Monte-Carlo 
method [2], the diffusion approximation [3,4], and the spherical harmonics discrete ordinate 
method [5-8]. 

In the last technique the phase function in the scattering integral is represented by the spherical 
harmonics and the integral is replaced by a quadrature sum. Spatial grids are introduced and 
obtained partial differential equations are approximated by the system of linear algebraic 
equations. To solve it, the successive-orders-of-scattering (SOS) approach is applied. The SOS 
approach has already been successfully applied in a number of cloud research topics, including the 
study of the transmission of thermal infrared radiation for a target—detector system [4] and the 
investigation of the effects of cloud geometry on the transmission of sunlight [9]. 

Each partial differential equation is integrated along its characteristic throughout whole 
calculation region in the framework of the well-known Evans’s algorithm [7]. Similar methods 
were developed for radiation transfer calculations in nuclear reactors shielding [10,11] and for the 
medical physics problems [12]. These methods have some advantages and some deficiencies. In 
particular, some of them may be not conservative, e.g. they do not conserve the number of 
particles (e.g., photons, electrons, neutrons) in the transport problem at hand. This defect can lead 
to significant errors in the solution obtained. 

Other spherical harmonics discrete ordinate methods use the local approximation for partial 
differential equations. Such methods were widely used in various neutron and photon transport 
problems in the last 50 years [3,13]. They are conservative and economic since these methods use 
very simple formulas and do not apply complicated logic (all spatial meshes are calculated 
successively). Namely, these methods have been incorporated in the RADUGA 3-D solver used in 
this paper [13,14]. They can be applied to a number of atmospheric radiative problems owing to 
two important features. Firstly, a special scheme for the partial differential equations solution, 
which obeys maximum principle, is used. Thereby, a grid solution conserves the exact solution 
positivity and all extrema dispositions are not corrupted by non-physical oscillations. Secondly, 
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this code can deal with multiprocessor computers and, therefore, can be used to calculate photon 
transport in large spatial regions. 

The RADUGA code can be applied to a number of atmospheric radiative problems, including 
studies of cloud and bright surface adjacency effects [15], cloud top topography/horizontal 
inhomogeneity effects [16,17], and biases in the cloud optical thickness t retrieval from space- 
borne platforms [18,19], to mention a few. 

In this paper the influence of a neighbouring cloud on the aerosol reflection function at the 
nadir observation geometry by means of the RADUGA code is studied. This is of importance for 
aerosol remote sensing (e.g., for a cloud screening procedures in satellite aerosol remote sensing 
algorithms [20]). 

The paper is structured as follows. In the next section of the paper we pose a physical problem. 
The third section is devoted to the introduction of the 3-D radiative transport equation. A brief 
discussion of the numerical procedure is also given there. A detailed description of the newly 
developed RADUGA code used in calculations is given in separate publications [13,21]. The 
numerical results obtained for the simple typical aerosol-cloud model are discussed in Sections 
4-6 of this work. 


2. The geometry of the problem 


The geometry of the problem is presented in Fig. 1, where a rectangular coordinate system xyz 
is introduced. Solar light is approximated by the monodirectional source having intensity Foô(u — 
M)6(o — ®). Here M = cos O, u = cos 0 and the pair (©, ®) gives the direction of solar light 
propagation in the spherical coordinate system defined by the axis z and angles (6@,@). The 
azimuth @ is counted with respect to the positive direction of the axis x. In this paper only results 
for ® = 0 and z will be reported. It means that solar light enters atmosphere from the direction of 


Satellite 


The line of 
visualization 


Fig. 1. The general geometry of the problem. 
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the positive values of x (then = r) or from the direction of the negative values of x (then ® = 0). 
The intensity of multiply scattered light is calculated along the axis x in the zenith direction as 
shown in Fig. | (see the line of visualisation in Fig. 1). 

We divide the terrestrial atmosphere in two equal semi-spaces separated by a local vertical 
plane. One part is filled by a cloudy medium and another one by the atmospheric aerosol. The 
processes of molecular scattering and absorption are neglected. Scattering media are assumed to 
be homogeneous and infinite in planes z = constant. We assume that there are no light scattering 
particles at z>4km. All downward propagated photons, which reach the plane z = 4, are 
assumed to be absorbed. Therefore, the contribution of the ground albedo is neglected. 

Droplets in a cloud are characterised by the Cloud C.1 particle size distribution [22] with the 
effective radius equal to 6um (see Appendix). The single scattering diagram for an elementary 
volume of a cloudy medium is calculated at the wavelength 412 nm using the Mie theory [23]. The 
phase function in the aerosol medium is represented by the Henyey—Greenstein formula [24] (see 
Appendix A). The asymmetry parameter of the cloud phase function g is equal to 0.85. The value 
of g for the aerosol phase function is equal to 0.7. The optical thickness of cloudy and aerosol 
portions of the scene are varied as specified below (see Section 4). Also we have studied the 
variation of the reflected light as observed from a satellite for a nadir observation geometry as the 
function of the solar angle ©. 


3. Theory 


The radiative transfer equation is used to solve the problem. This equation has the following 
form for the case studied [13]: 


Ol(x, y, Z, 8, p) Ol (x, y, Z, 8, p) 
+ 
Ox Oy 


1 T f 2n ; , 
E Gscal x, y, Z) f sin 0’ d0 f I(x, y,z,9', o')p(x, y, z, 18, 9, 0, o) dg’ 
0 0 


I(x, y, Z, 8, p) 
Oz 


č + B + Cext(X, Y, Z)I(x, y, Z, 0, p) 


1 
T dn Osca(X, y., z)Fop(x, J; Z, x(0, Ọ, O, ®)) exp(—?), (1) 


where osca and Gex, are scattering and extinction coefficients, which are assumed to be identical in 
this study (possible light absorption both by droplets and aerosol particles is neglected). 
The function /(x,y,z,0,@) is the diffused light intensity at the point (x,y,z) propagated 
in the direction (0, p), see Fig. 1. Also we have: € = sin @cos g, n = sin Osin ọ, P = cos 0. The 
phase function p(x, y, z, X) is defined in Appendix A. Scalar product x is defined by the following 
relation: 


x(0, p, 8,0’) = cos O cos 6’ + sin Osin 6 coslo — ^). (2) 
ne value of t in Eq. (1) is the optical path between two points defined by radius-vectors 79 and 


Ei i Gext(7o + 67) df, where d = |F — ro|. The vector Fo defines the cross point of the light 
sean with the boundary of scattering medium under study. 
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We assume that there is no diffused light entering the medium. Therefore, boundary conditions 
have the following form: 


1(?,Q) =0 at OAP) <0 for all F(x, y, Z) € Grnites 


1(7,Q) = 1(7",Q) at Q(#) <0 for all Xx, y, z) € Gintinite- 


Here Gfinite is the finite part of the medium boundary (on z), Ginfinite is the infinite one (on x and y), 
r“ is an inner point for the periodic boundary conditions. 

We will assume that the line LM (see Fig. 1) coincides with the visualisation line. It means that 
the solution becomes invariant in respect to the coordinate y. Then we can drop the dependence 
on y in Eq. (1) and arrive to the following simplified 2-D transport equation (see Fig. 2): 


I(x, z, 0, p) +B I(x, z, 0, p) 


+ Cext(X, Z)I(x, Z, 0, p) = FI, (3) 
Ox Oz 


č 


where 


T 2m 
FI = - Gacy XZ) f sin 0’ dé’ f I(x, z, 0’, p’)p(x, z, x0, 0, 0, PN) dq’ 
T 0 0 


1 
+ ge z)Fop(x, Z, (0, P, O, ®)) exp(— ô). (4) 


Although note that the RADUGA can be applied to arbitrary-shaped scattering media. 

Eq. (3) is solved using the method of successive orders of scattering. Namely, at first we neglect 
the integral term in Eq. (4) and calculate the diffused intensity /(x, z, 0, p) from the solution of the 
partial differential equation. Then the obtained diffused intensity is substituted in the scattering 
integral in Eq. (4) and the next approximation for I(x, z, 0, p) is found from the solution of the 
partial differential equation (3). The algorithm is stopped when the convergence is reached. 

Therefore, the problem at hand is reduced to the solution of the following transport equation: 


I(x, z, 0, p) I(x, z, 0, p) 
Ox Hi Oz 


where I(x,z,0, p) is the diffused light intensity at the point r(x,z) propagated in the direction 
(6,~), I(x, z,0, p) = FI(x,z,8,~), and I(x,z,8,@) is the function obtained from the previous 


iteration as described above. 
rhe gf y fOlY K Line of 
OTA D| C visualization 
l 


¢ + Cex X) (x, Z, 0, p) = I(x, Z, 0, p), (5) 


Aerosol q—> ~ Cloud 
I I 
l I 
7 O i , 
0 x, km 
z, kmy 


Fig. 2. The geometry of the 2-D problem. 
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We introduce an angular quadrature and replace functions /(x, z, 0, ¢), I(x, z, 0, p) and 
(x, z,9,@) by their values in quadrature nodes. The integral 3(x,z,0, ø) is calculated using 
following standard steps: 


@ the expansion of the function I(x, z, 6, ọ) in terms of spherical harmonics, 
@ the representation of the integral X(x, z, 0, p) by spherical harmonics, 
è the calculation of 3(x, z, 0, @) values in quadrature nodes. 


The details will be given elsewhere [13] and will be not considered here. Variables 0 and @ are 
omitted for simplicity. 

A standard grid method to approximate the partial differential equation of the first 
order (5) is used. In particular, grids with the respect to spatial variables x and z are 
introduced: 


X12 < ttt SX 2S SX K41/2, Z1/2< tt SZ 1/2< +++ KZL+1/2- 


A single two-dimensional cell (k,/) has following dimensions: [x;~1/2, Xx41/2] x [2-1/2 2141/21: 
Correspondingly, its size is [Ax,] x [Az;], where Ax, = Xk+1/2 — Xk-1/2 and AZp = Zb4+1/2 — Zb—1/2- 
Also the integral operator 


p 1 Xk+1/2 Z1+1/2 
Ax; = —— f dx f dz 
AxXAy] Xk-1/2 Z]-1/2 


is applied to both parts of Eq. (5). Then it follows: 


k,l ~ 
ElI k+i/21 — Lk-1/21)/Axp + BU 2 — Let1/2)/Az1 + Cele = Sk (6) 
where 
1 Xk+1/2 Zl+1/2 1 Xk+1/2 Zl+1/2 
Iki = zx] dx f dz I(x, z), Ik, = / dx | dz Sx, Z) 
Ax, Az; Xk-1/2 2-1/2 Ax. Az] Xk—1/2 Z1-1/2 
are the average values of the intensity and the source function, respectively, over a given cell and 
1 Z141/2 1 Xk+1/2 
Iks = — dz I(xk+1/2;Z),  Ikiz1/2 =- dx I(x, Z+ 
ke1/2d = Az L. Z1(Xk41/2.2),  Iki41/2 Dr Jaan (x, 2741/2) 


are correspondent average values of the intensity on boundaries of the cell (k, 1). Fulfilment 
of Eq. (6) guarantees that the presented scheme is a conservative one. 

Intensities 7,_/2, and Ik -1/2 are known either from boundary conditions or from the result of 
the calculation for the previous cell. So we need to determine only values of Ix41/27, Zki+1/2 and 
Iķı. It is not possible to evaluate three parameters from a single equation (6). So we need to 
introduce two approximate relations among these three unknown parameters. They are given as 
follows: 


Iki = (l = ty) Ter s(Q20 + Vxk H k-s8/21 


Iki = (1 — vk) ki+5(p)/2 + Vk, k,I-s8)/2> 
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where s(¢) = sign(€), s(p) = sign(P), vx kı € [0, 1) and vz; € [0, 1) are weight parameters. We use 
values of weight parameters as introduced in [25,26]: 


_ (hy ii /hz kA + hxk) 
2+ 2hyx k + ae 


tego = 1/(2 + hxk), Vak at hyk I <hz x, 


(he pt /hy kA + hz) 
c= 5 
2 + 2h-k1 + IK) 


» Vaki = 1/(2 + h-ki) at hye Shek. 


Optical steps xx 1 and h,,, are defined as 


k,l k,l 
hki = Cot Axe cl. haki = onze) |p. 


This scheme defines outside fluxes Ik+s/2,1 Tki+s(p)/2 Via entering fluxes Ik-s6/2,1; Ik 1-s(8)/2 1N a 
physically correct manner in any cell of any grid [25,26]. It permits to obtain discontinuous 
solutions and ones with great gradients with a high accuracy. This completes a brief description of 
our technique. Further details are given by Bass et al. [13,14,25,26]. 


4. The validation of the RADUGA code 
4.1. The comparison with the 1-D transport code 


Under conditions specified above, the considered problem is reduced to the 2-D problem 
presented in Fig. 2. Both aerosol and cloud are homogeneous along z-direction on the height 
interval (0 km, 4km). They are contained in rectangular boxes. The line OD (see Fig. 2) separates 
aerosol medium from a cloud. Boxes are of an infinite length in the direction perpendicular to the 
plane of the figure of theand along axis x. 

We will study the upwelling light field in the zenith direction at line AC. Clearly, the intensity of 
reflected field must depend on the coordinate x. This corresponds to the case of an orbiting 
satellite or aircraft observing aerosol-cloud system from above. 

The reflection function R is defined as the ratio of the light intensity Z reflected from a given 
medium to that reflected from an absolutely white Lambertian surface /,. 


R=I/I,. (7) 


It is easy to show that J, = cos © Fo/z, where Fo is the density of incident flux on the unit area 
perpendicular to the beam and @ is the incident angle [24]. 

The largest gradients of the function R(x) are expected in the area closest to the cloud boundary 
OD. Because both an aerosol medium and a cloud are extended to infinity along axis x, this 
function far from boundaries must be equal to the value, which can be obtained from the 1-D 
radiative transfer equation. 

To validate our results we calculate the dependence R(x) at large distances from the boundary 
OD using the RADUGA code and compare results obtained with values derived using an 
independent 1-D code based on the discrete ordinates method as described in [27]. These results 
are affirmed by the asymptotic theory for optically thick slabs [28,29] in the case of a cloud. The 
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Table 1 
The comparison of the 3-D RADUGA code with 1-D code [27] results at large distances from the aerosol-cloud 
boundary 


Aerosol (x + —oo) Cloud (x + oo) 
(8) 1-D 3-D 1-D 3-D 
50 0.1201 0.1201 0.6921 0.6925 
60 0.1494 0.1494 0.6493 0.6493 
70 0.1822 0.1823 0.5666 0.5638 
80 0.1977 0.1976 0.4476 0.4474 


comparison is shown in Table | assuming that the aerosol optical thickness t* is equal to 1.2 and 
the cloud optical thickness t° is equal to 30. The solar angle is varied in the range 50—80°. 

It follows from the analysis of Table 1 that the difference of results obtained from different 
codes is below 0.5% in the aerosol medium. It is smaller than 1% in the cloud. Therefore, we 
conclude that RADUGA is capable to reproduce 1-D code results with a high accuracy. 


4.2. The comparison of results obtained from the RADUGA code with the Monte-Carlo calculations 


We also performed extensive comparisons with the Monte-Carlo calculations based on the 
MYSTIC (the Monte Carlo code for the physically correct tracing of photons in cloudy 
atmospheres, see [30,31]). The MYSTIC is a forward Monte Carlo code which traces photons on 
their individual paths through the atmosphere, similar to what is described in [32]. Radiances are 
calculated using a local estimate technique (e.g., [2,33]). In this configuration, the MYSTIC has 
been successfully validated in the intercomparison of 3D radiation codes (see http:// 
climate.gsfc.nasa.gov/I3RC). The MYSTIC is operated within the libRadtran package (see 
http://www.libradtran.org), which prepares the optical properties of the atmosphere, to be used in 
the model. For the application in this paper, the atmosphere was consisted of one layer, with 
aerosol and cloud properties as specified in this paper. A model domain of 80 km in x was used. A 
large domain size is important since the MYSTIC uses periodic boundary conditions. The model 
resolution was set to 0.1 km; the MYSTIC results are, therefore, averages over 0.1 km bins. 

Some results of comparisons are given in Figs. 3a and b (at O = 60° and azimuths # = 0 
and z). We find that differences are below 1%. Therefore, the RADUGA provides very 
accurate results as far as calculations of the light reflected by an aerosol-cloudy medium are of 
concern. 


5. Main properties of computed reflection functions 


Some physical dependencies are clearly seen in Figs. 3a and b. For instance, it follows from 
Fig. 3b that there is a shadow near the cloud border for the illumination from the cloud side 
(@ = z). Also we have a brightening effect in Fig. 3a due to the cloud side illumination effects 
(@=0). These two effects (shadowing and brightening) are primarily due to the direct light 
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Fig. 3. The dependence R(x) calculated using the Monte-Carlo code MYSTIC and the RADUGA code for the case 
shown in Fig. 2 at the illumination from the cloud side (a) and the aerosol side (b) for the solar angle 60°. The aerosol 
optical thickness t* is equal to 1.2 and the cloud optical thickness t° is equal to 30. The phase functions are specified in 
the Appendix. The absorption of light and the contribution from the surface reflection is neglected. The cloud boundary 
is placed at x = 40 km. 


interaction with a scattering medium. They lead to roughening effects in 2-D—3-D transfer 
problems. We also observe (see Fig. 3b) the decrease of the reflection near the border of the cloud 
(inside the cloud) as compared to 1-D case. This is due to photon leaking in the area with smaller 
extinction coefficient in the direction of the part of the scene. The increase in the aerosol reflection 
function in Fig. 3b close to the cloud is due to channelling of photons from a cloud to the aerosol 
side. These two effects (photons channelling and leaking) lead to the smoothening of the radiative 
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field. Four effects considered here exist not only for a simple case studied in this paper but also for 
broken cloud systems [34,35]. 

The dependence R(x) can be easily parameterised taking into consideration these four 
fundamental phenomena. The parameterisation of 3-D effects is of a great importance for satellite 
remote sensing because calculations presented here are computationally expensive and cannot be 
included in the operational aerosol/cloud retrieval algorithms [29]. Also these effects cannot be 
neglected. In particular, if an orbiting optical instrument observes an area correspondent to 
brightening/shadowing effects, then large biases in retrieved cloud/aerosol parameters are 
expected. This also points to the necessity of the development of the simultaneous aerosol—cloud 
(and surface) retrieval algorithm. The complex system should be considered in the retrieval as a 
whole [36]. 


6. The dependence of the aerosol reflection function on the solar angle 
6.1. The illumination from the cloud side 


1. The function R(x) for the aerosol-cloudy medium specified above is shown in Fig. 4 for the 
case T? = 0.25, t° = 30 . This dependence was obtained using the RADUGA code for solar angles 
@ in the range 10—-80°, = xz, and the nadir observation. The characteristic feature of functions 
given in Fig. 4 is the presence of the brightness minimum in the aerosol area close to the cloud 
boundary. Because the aerosol is optically thin, the aerosol reflection function chiefly depends on 
the singly scattered component R;(x). 


Reflection function 


Fig. 4. The same as in Fig. 3 but for the aerosol optical thickness 0.25 and various solar angles (10-80°) at the 
illumination from the cloud side. 
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The equation for the function R; can be derived from the relations (2)-(4), (7) at 0 = 2, any ọ, 
and I(x, z, 0’, o) = 0. Namely, it follows y = — cos @ and 
OR, Oext(X, z) 


= gr TIR, 2) = EO P(x, z,—cos O)e"™),  —co<x<00, O<z</. (8a) 


Here / = 4km is the cloud top height. The optical path ¢(x,z) is determined by the following 
formula: 


o2,42/ cos © as z< — xctgO, x<0, 
t(x,z) = l 9542z/ cos O — [ot — of,,)x/sinO@ asz>—xctgO, x<0, 
a Z/ COS O as x>0, 


where oĉ and of, are the extinction coefficients for the aerosol and the cloud, respectively. The 
boundary condition corresponds to the absence of the light propagated upwards at the border 
z= l: 


R(x, D = 0. (8b) 


The problem (8) has a unique solution. Namely, it follows: 


_ P&Z- cos 8) P , 
Ry (3,2) = PO ala) | PO DE- D expl O) A, 
Substituting z = 0, it follows: 
Ri(x) — u(O){1 _ e/4*(@)/ cos ©] as X<X,, (9) 


R(x) — u(@){1 + qo! (Oje ONO) = d*(O)v!(O)e)/ 8 O e0) as x, <x<0, 


R(x) = w(O)[1 — e/s e] as x>0. 


Here 
x, = —ItgO, (10) 
s(x,0)=x/sin ©, u”(@)= p”(— cos O)/[4(1+ cos O)], q = dha — Cap (la) 
u(O) = oba cos O + 0go d”(O) = oh (1 + cos O), (11b) 


where m = a,c and p* and p° are the aerosol and cloud phase functions, correspondingly. 
The function R(x) is constant as x>0 and as x<x;. In interval [x,,0) it is smooth and its 
derivative is given by following equation: 


(d/dx)Ri(x) = qê (OWT (Od (cos Oje” OOP] — AOE cos 1 sin O. 


Since for the considered problem the value q is negative, then the component R(x) monotonically 
decreases in this interval. The component R(x) and the corresponding full solution R(x) at the 
angle O = 60° are plotted in Fig. 5. 

At the media boundary (point x = 0) the function R(x) has a discontinuity and the value Rı(0) 
is assumed to be indefinite [37]. This discontinuity of the singly scattered component Rj(x) is due 
to the discontinuous optical properties of the aerosol-cloud medium under study. It is diminished 
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Reflection function 


Fig. 5. Reflection function R(x) and its singly scattered component for solar angle 60° at the illumination from the 
cloud side. 


by scattering acts and forms minimum of the full solution R(x), look in Fig. 5. This minimum is 
the manifestation of the cloud shadowing effect. 

2. An important point in the aerosol remote sensing from space is the estimation of the distance 
from the cloud, where the aerosol optical thickness can be retrieved using the 1-D transport 
theory. To answer on this question we present data given in Fig. 4 in yet another form. 

Figs. 6a and b show the dependence of the relative difference 6 (in percent) between the 
reflection function R(x) obtained from the 3-D code as compared with 1-D calculations R* for the 
aerosol layer far from the cloud boundary. Namely, we define 


O(x) = 1 — R(x)/R*. 


Clearly, it follows as x + —oo: 6(x) > 0. Circles in Fig. 6b show points, where x is equal to xs. 
The value of x, gives the exact coordinate of the cloud shadow boundary for singly scattered 
component R(x), see Eq. (10) and Fig. 5. Owing to scattering processes for the full function R(x) 
the value x, is the estimate from below for the cloud shadow boundary coordinate. 

We see that the functional dependence 6(x) changes its slope at x< xs. If we tolerate the 15% 
error in the value of R, then we can state that 1-D theory can be applied at x< x, for the case 
studied here. Note that for the quantitative aerosol remote sensing values of R should be correct 
at least within 5%. 

It means in particular, that the influence of cloud is of importance at distances 4.8, 6.9, 11.0, 
22.7km for angles equal to 50°, 60°, 70°, 80°, respectively, for the case studied here. Clearly, no 
aerosol optical thickness retrievals are possible at the asymptotic case O = 90° (x, —> —oo). It is 
interesting to see that the value of ô is positive for large solar angles. This means that 2-D-3-D 
effects decrease the reflection function as compared to a simple 1-D theory in the case under 
study. Surprisingly, for sun close to the zenith the light enhancement reflectance in the clear sky 
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Fig. 6. The dependence ô(x) in the clear-sky part of the scene for solar angles 10—40° (a) and 50-80° (b) degrees 
correspondent to the case shown in Fig. 4. 


region is possible (see Fig. 6a). For solar angles smaller than 30° and for a 10% threshold the 
influence of cloud is of importance up to the distance 10 km. This threshold is somewhat smaller 
(approximately, 5km) for solar angles 40—50° but substantially increases for large solar angles (see 
Fig. 6b). 

We also can plot the function ô(x) for a cloudy part of the scene. The dependence of the relative 
difference ô close to the cloud boundary is monotonous in this case (see, e.g., Fig. 7). This is 
primarily due to the absence of a shadowing effect for a cloud. Also we note that ô(x) is only 
weakly dependent on the solar angle. In particular, data presented in Fig. 7 can be approximated 
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Fig. 7. The same as in Fig. 6 but for the cloudy part of the scene. 


by the following curve: 
0(x) = A exp(—x/X0), (12) 


where we introduced the correlation length xo. It gives the distance, where the difference 
attenuates in e times. We find from data given in Fig. 7 that A = 54.4 and x9 = 1.4km. The value 
of xo should correlate with the transport length ly =1/(1 — g), where 1 = //t is the photon free 
path length in the cloud. It follows for the case we consider: 1 ~ 0.133 km7!, g ~ 0.85, and, 
therefore, ly ~ 0.889 km. So we have: xp © 1.6l;;. The distance, where the difference 6 is below 
10% is approximately equal to 3l in the case studied. 

3. Yet another effect is seen in the reflection function R(x). Its constant value in the aerosol area 
far from the media boundary is more for the low sun (close to the horizon) as compared to the 
case of a high sun. Similar assertion is true for the singly scattered component. This follows from 
the expression for the derivative of the function R(x; @) with respect to the parameter © 


sin O 
p*(— cos O) 
+ Rj(x)[(0/0 cos O)(p*(— cos @)) + 4u*(O)]} as x<xs, 


the positivity of Henyey—Greenstein phase function derivative and an inequity u*(@) — R(x) 20, 
see Eq. (9). 

The opposite effect is observed for the reflection function in the cloud: its value is smaller for a 
high sun position as compared to the case of a low sun. This effect cannot be explained by singly 
scattered component analysis because the cloud is optically thick and the function R(x) is formed 
by many scattering acts. To clarify this feature, we note that if a photon entered an optically thick 
non-absorbing cloud through a point of its top border placed far from the other borders, then this 
photon will leave the cloud through the top boundary after many collision events. The more 


(d/d@)R\(x; @) = {p*(— cos @)[u*(@) — Ri(x)]/o°,, /cos” O 
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collisions the photon undergoes the more its outcoming direction differs from the incident one. As 
a consequence, it follows that the more an incident direction differs from the axis z the smaller 
upward reflection function value in the direction perpendicular to the cloud border. 


6.2. The illumination from the aerosol side 


The results of calculations for the illumination from the aerosol side (® = 0) are given in Fig. 8. 
The sequence of curves in the aerosol/cloud side outside the transition region is similar to that 
given in Fig. 4 and has the same explanation. A new effect is seen however. Namely, it follows that 
the cloud reflection function increases close to the cloud boundary. This is due to the illumination 
of a cloud from a side. Therefore, the interaction of direct light with a cloud may lead not only to 
a familiar shadowing as in Fig. 4 at x ~ 0 but also to the brightening effect depending on the 
cloud position with respect to the incident light direction. 

Such effect may be explained by the singly scattered component consideration. This function 
for this case is defined by relations 


R(x) = u(O){1 = e /4(0)/ cos 0) as X> Xy, 
RI) = (O1 — e’ Oos 0) as <0, 
Rı(x) r wW(OH1 — qu (O)e TOO) 


_ d°(O)v!(O)e)/ 08 Peas) as 0<x<x,. 


Here x, = Ltg O and other parameters are determined by the formulas (11). 
In the case of illumination of cloud from the aerosol side, the dependency of the aerosol 
reflectance function on x is a monotonous one (see Figs. 9a and b). Generally, the aerosol 


Reflection function 


Fig. 8. The same as in Fig. 4 but the illumination from the aerosol side. 


420 O.V. Nikolaeva et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 94 (2005) 405—424 


-100 


—— © = 30° and © = 40° 


N 
[s] 
ls} 


8 


Difference 8(x) (%) 


(a) 


Difference (x) (%) 


8 8s 8 & 8 o 


(b) 


Fig. 9. The same as in Fig. 5 but for the illumination from the aerosol side. 


brightness increases in the vicinity of the cloud. This increased brightness of the aerosol may lead 
to the positive bias in the retrieved aerosol optical thickness as derived from data obtained by 
space born radiometers and spectrometers. Therefore, aerosol retrievals should be performed at 
least at the distance 10-15km from the cloud boundary (see Figs. 9a and b) depending on the 
illumination conditions. 

Unlike the case of the illumination from the aerosol side, the cloud reflection function has a 
maximum, which is due to the brightening effect. The brightening effect complicates the cloud 
optical thickness retrieval from space. In particular, to avoid this effect only pixels positioned at 
x > 2-8 km (depending on the solar angle, see Fig. 10) should be used in the retrieval procedure. 
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Fig. 10. The same as in Fig. 6 but for the illumination from the aerosol side. 


Therefore, we see that 2-D—3-D effects influence considerably not only the aerosol reflectance 
but also the same is true for a cloudy part of the scene under study. 


7. Conclusions 


The developed 3-D radiative transport code RADUGA has been tested and applied to studies 
of the influence of cloud boundaries on the clear sky reflectance function. In the considered 
boundary problem for the transport equation we have investigated the phenomenon of the inner 
boundary layer with transient radiation regimes. The reflectance properties of this boundary layer 
are studied in details. We found that the aerosol reflectance is sensitive to neighbouring clouds for 
distances up to 25km for realistic observation scenarios, depending on the solar angle. In this 
paper the cloud top height was fixed. Clearly, higher clouds can influence the aerosol zenith 
reflectance even at larger distances. Therefore, special caution should be taken in clear sky aerosol 
optical thickness retrieval from space even if it seems that clouds are far away from a scene 
studied. The same applies to the cloud optical thickness retrieval near cloud boundaries (e.g., see 
Fig. 10). 

We conclude that 3-D effects substantially reduce the number of pixels, which can be used in 
the retrievals based on the 1-D radiative transfer theory. To eliminate this difficulty, retrievals 
should be based on the 3-D radiative transfer theory. However, this is difficult to achieve on the 
operational basis due to the limitations imposed by modern computers. Therefore, there is an 
urgent need for the parameterisation of 3-D effects in cloudy atmospheres. Such parameterisa- 
tions can be easily performed for regions where radiative smoothing takes place. Corresponding 
dependencies 6(x) can be approximated by rather simple exponential functions (see, e.g., Eq. (12)). 
Parameterizations become a little bit more complex for shadowing/brightening areas due to a 
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non-monotonous behaviour of studied functions there. Also the dependence on the solar angle is 
much more pronounced then (compare, e.g., Figs. 6 and 9). 

We note that other radiative characteristics (e.g., the cloud transmittance and fluxes) can be 
also studied in the framework of the code developed. 
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Appendix A. Phase functions 


The phase function p(z) gives the conditional probability of light scattering by an elementary 
volume of a scattering medium in dependence of the cosine x of the scattering angle (the value 
y= 1 corresponds to the forward scattering). In the transport calculations this function is 
represented usually by the expansion in terms of Legendre polynomials P,,(y) 


N 1 1 
1G) = S°On+ DP m= POPA (A.1) 
n=0 


where, in particular, Po) = 1, P10) = x, P20) = (37° — 1)/2. The constant g = g; is called the 
asymmetry parameter. 

The phase function of light scattering by spherical water droplets in clouds can be calculated 
using the Mie theory. The dependence of the phase function on the cosine of the scattering angle x 
for the Cloud Cl model [20] at the wavelength 412nm is shown in Fig. 11 (left panel). Its 
asymmetry parameter is close to 0.85. We use this phase function to represent light scattering by 
an elementary volume of a cloudy medium in this paper. 
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Fig. 11. The aerosol and cloud phase functions. 
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Equality (A.1) for the cloud phase function has a high accuracy at N>350. The number of 
terms N influences the speed of calculations considerably. Therefore, we applied the delta-M 
method [38] with the aim of the reduction of number of terms in the correspondent expansions. 
Basically, in this method the phase function is presented as a sum of a singular part (represented 
by a delta function) and a regular part, which can be represented by the combination of much 
smaller number M of Legendre functions (M < N). In this paper we have used M = 28 (see Fig. 
11). It follows that the resulted phase function (dashed line) closely corresponds to the initial Mie 
phase function (solid line). Note that the radiative transport in thick clouds (due to strong 
multiple scattering) is determined by a few first moments of the phase function (e.g., the 
asymmetry parameter g, see above). It is insured that 27 moments of the exact and approximated 
phase functions coincide in our case. We believe that such an approach does not influence the 
reflection function calculations in practical terms, which is also confirmed by comparison with 
other codes. 

The phase function of atmospheric aerosol cannot be calculated using the Mie theory in the 
general case. This is due to the non-spherical shape of a great portion of aerosol particles (e.g., dust). 
So it is assumed that the aerosol phase function can be presented by Henyey—Greenstein formula: 


pO) = (1 — 9)/0. +P — 29x)”, 


here g, = g”. We choose g = 0.7 (such value is a typical case for atmospheric aerosol) and assume in 
practical calculations that N = 23. The accuracy of this approximation is shown in Fig. 11 (right 
panel). The error of approximation is negligible in practical terms. 

Although this simplified aerosol phase function may be unrealistic as far as real atmospheric 
aerosol is of concern, it does not influence the main result of this work, namely the cloud adjacency 
effects studies. Also other researches can easily incorporate this phase function in their codes, 
which can facilitate the comparison of our results with those of others. We also neglect the vertical 
variability of the phase function as not an essential issue for the problem studied in this paper. 
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